
#### Preliminaries ####

R.version$version.string

# clear workspace
rm(list=ls())


# load packages
library(here)
library(tidyverse)
library(rworldmap) # maps
library(magrittr)
library(reshape2)
library(ggpubr)
library(leaflet)
library(plotly)
library(gapminder)
library(countrycode)
library(ggrepel)
library(ggthemes)
#library(ERT)
library(modelsummary)
library(zoo)
library(brglm)
#devtools::install_github("vdeminstitute/ERT", dependencies=TRUE)
library(ERT)
library(DIDmultiplegt) 


#---------------------------------------------------#
#---- Load Function for getting Episodes of X -----#
#---------------------------------------------------#

# load the relevant functions into the environment

source("R/get_episode_wo_CI_function_four_item_afi.R")

#### Data merging ####

vdem_full_v13 <- readRDS("data/vdem/V-Dem-CY-Full+Others-v13.rds")

academic_freedom_index_data_four_item <- readRDS(file.path("results/results/cy", "bayes.academic_freedom_index_cy.rds"))

academic_freedom_index_data_four_item <- academic_freedom_index_data_four_item %>%
  rename(v2xca_academ_four = mean, 
         v2xca_academ_sd_four = sd, 
         v2xca_academ_codelow_four = codelow68, 
         v2xca_academ_codehigh_four = codehigh68) %>%
  dplyr::select(country_text_id, year, ends_with("four"))

vdem_full_v13 <- vdem_full_v13 %>%
  left_join(academic_freedom_index_data_four_item, by = c("country_text_id", "year"))

vdem <- vdem_full_v13

episode_data_wo_CI <- get_episode_wo_CI_function_four_item_afi.R()

rm(vdem)

ert_episode_data <- read.csv("data/ert.csv")

#### prepare merging ####

ert_episode_data <- ert_episode_data %>%
  dplyr::select(-c(country_name, country_id, v2x_regime, starts_with("v2x_polyarchy")))

episode_data_wo_CI <- episode_data_wo_CI %>%
  ungroup() %>%
  dplyr::select(- c(v2xca_academ_four, country_name, country_id))

#### merge data ####

vdem_subset_onset_sample <- vdem_full_v13 %>%
  left_join(episode_data_wo_CI, by = c("country_text_id", "year")) %>%
  left_join(ert_episode_data, by = c("country_text_id", "year"))

# Create a tidier for "multiplegt" objects
tidy.did_multiplegt = function(x, level = 0.95) {
  ests = x[grepl("^placebo_|^effect|^dynamic_", names(x))]
  ret = data.frame(
    term      = names(ests),
    estimate  = as.numeric(ests),
    std.error = as.numeric(x[grepl("^se_placebo|^se_effect|^se_dynamic", names(x))]),
    N         = as.numeric(x[grepl("^N_placebo|^N_effect|^N_dynamic", names(x))])
  ) |>
    # For CIs we'll assume standard normal distribution
    within({
      conf.low  = estimate - std.error*(qnorm(1-(1-level)/2))
      conf.high = estimate + std.error*(qnorm(1-(1-level)/2))
    })
  return(ret)
}
theme_set(theme_minimal()) # Optional


############################################################################################
#### Descriptive Statistics                                                             ####
############################################################################################

vdem_descriptives <- vdem_subset_onset_sample %>%
  dplyr::select(country_name, year, v2xca_academ_four, v2x_polyarchy, dem_ep, aut_ep, 
                increase_episode, decline_episode, 
                e_gdppc,e_pop, v2x_jucon, v2xlg_legcon) %>%
  drop_na(v2xca_academ_four)

vdem_descriptives$year <- as.integer(vdem_descriptives$year )
vdem_descriptives$country_name <- as.factor(vdem_descriptives$country_name)


datasummary_skim(vdem_descriptives, 
                 title = "Descriptive Statisitcs Sample", 
                 output = 'tables/Table_F2.tex', 
                 fmt="%.3f",
                 histogram=FALSE)

summary(vdem_descriptives)

############################################################################################
#### Figure F1: Density Plot Both AFI indices                                            ####
############################################################################################


ggplot(vdem_subset_onset_sample) +
  geom_density(aes(v2xca_academ), color = "red") +
  geom_density(aes(v2xca_academ_four), color = "blue") +
  theme_bw() +
  xlab("Academic Freedom") +
  ylab("Density") 

ggsave("figures/Density_AF.pdf")



############################################################################################
#### Figure 6: Analyzing Democratization's Effect on AF                                 ####
############################################################################################

m1 = did_multiplegt(
  vdem_subset_onset_sample, 
  Y = 'v2xca_academ_four', 'country_id', 'year', 'dem_ep', 
  dynamic   = 10,                  # no. of post-treatment periods
  placebo   = 10,                  # no. of pre-treatment periods
  brep      = 20,                  # no. of bootstraps (required for SEs)
  cluster   = 'country_id',        # variable to cluster SEs on
  parallel  = FALSE                 # run the bootstraps in parallel
)


tibble(estimate = m1$effect)


tidy_m1 = tidy.did_multiplegt(m1)
tidy_m1

tidy_m1 |>
  within({
    term = gsub("^placebo_", "-", term)
    term = gsub("^effect", "0", term)
    term = gsub("^dynamic_", "", term)
    term = as.integer(term)
  }) |>
  ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange() +
  geom_hline(yintercept = 0, linetype = 'dashed') +
  labs(
    x = "Time to treatment (Democratization Episode)", y = "Effect size", 
    title = "", 
    subtitle = ""
  ) 

ggsave("figures/Figure_1_Demo_Effect_on_AF.pdf", height = 10, width = 20, units = "cm", dpi = 900)
ggsave("figures/Figure_1_Demo_Effect_on_AF.png", height = 10, width = 20, units = "cm", dpi = 900)


############################################################################################
#### Figure F2: Analyzing Democratization's Effect on AF                                ####
############################################################################################

m1b = did_multiplegt(
  vdem_subset_onset_sample, 
  Y = 'v2xca_academ_four', 'country_id', 'year', 'dem_ep', 
  controls = c("e_gdppc","e_pop", "v2x_jucon", "v2xlg_legcon"),
  dynamic   = 5,                  # no. of post-treatment periods
  placebo   = 10,                  # no. of pre-treatment periods
  brep      = 20,                  # no. of bootstraps (required for SEs)
  cluster   = 'country_id',        # variable to cluster SEs on
  parallel  = FALSE                 # run the bootstraps in parallel
)


tidy_m1b = tidy.did_multiplegt(m1b)
tidy_m1b


tidy_m1b |>
  within({
    term = gsub("^placebo_", "-", term)
    term = gsub("^effect", "0", term)
    term = gsub("^dynamic_", "", term)
    term = as.integer(term)
  }) |>
  ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange() +
  geom_hline(yintercept = 0, linetype = 'dashed') +
  labs(
    x = "Time to treatment (Democratization Episode)", y = "Effect size", 
    title = "", 
    subtitle = ""
  )

ggsave("figures/Figure_1_Demo_Effect_on_AF_with_covariates.pdf", height = 10, width = 20, units = "cm", dpi = 900)
ggsave("figures/Figure_1_Demo_Effect_on_AF_with_covariates.png", height = 10, width = 20, units = "cm", dpi = 900)


############################################################################################
#### Figure 7: Analyzing Autocraitzation's Effect on AF                                 ####
############################################################################################

m2 = did_multiplegt(
  vdem_subset_onset_sample, 
  Y = 'v2xca_academ_four', 'country_id', 'year', 'aut_ep', 
  dynamic   = 10,                  # no. of post-treatment periods
  placebo   = 10,                  # no. of pre-treatment periods
  brep      = 20,                  # no. of bootstraps (required for SEs)
  cluster   = 'country_id',        # variable to cluster SEs on
  parallel  = FALSE                 # run the bootstraps in parallel
)


tidy_m2 = tidy.did_multiplegt(m2)
tidy_m2


tidy_m2 |>
  within({
    term = gsub("^placebo_", "-", term)
    term = gsub("^effect", "0", term)
    term = gsub("^dynamic_", "", term)
    term = as.integer(term)
  }) |>
  ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange() +
  geom_hline(yintercept = 0, linetype = 'dashed') +
  labs(
    x = "Time to treatment (Autocratization Episode)", y = "Effect size", 
    title = "", 
    subtitle = ""
  )

ggsave("figures/Figure_2_Auto_Effect_on_AF.pdf", height = 10, width = 20, units = "cm", dpi = 900)
ggsave("figures/Figure_2_Auto_Effect_on_AF.png", height = 10, width = 20, units = "cm", dpi = 900)

############################################################################################
#### Figure F3: Analyzing Autocraitzation's Effect on AF                                ####
############################################################################################

m2b = did_multiplegt(
  vdem_subset_onset_sample, 
  Y = 'v2xca_academ_four', 'country_id', 'year', 'aut_ep', 
  controls = c("e_gdppc","e_pop", "v2x_jucon", "v2xlg_legcon"),
  dynamic   = 5,                  # no. of post-treatment periods
  placebo   = 10,                  # no. of pre-treatment periods
  brep      = 20,                  # no. of bootstraps (required for SEs)
  cluster   = 'country_id',        # variable to cluster SEs on
  parallel  = FALSE                 # run the bootstraps in parallel
)


tidy_m2b = tidy.did_multiplegt(m2b)
tidy_m2b


tidy_m2b |>
  within({
    term = gsub("^placebo_", "-", term)
    term = gsub("^effect", "0", term)
    term = gsub("^dynamic_", "", term)
    term = as.integer(term)
  }) |>
  ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange() +
  geom_hline(yintercept = 0, linetype = 'dashed') +
  labs(
    x = "Time to treatment (Autocratization Episode)", y = "Effect size", 
    title = "", 
    subtitle = ""
  )

ggsave("figures/Figure_2_Auto_Effect_on_AF_with_covariates.pdf", height = 10, width = 20, units = "cm", dpi = 900)


############################################################################################
#### Figure 8: Analyzing Academic Freedom Growth's Effect on Democracy                  ####
############################################################################################

m3 = did_multiplegt(
  vdem_subset_onset_sample, 
  Y = 'v2x_polyarchy', 'country_id', 'year', 'increase_episode', 
  dynamic   = 10,                  # no. of post-treatment periods
  placebo   = 10,                  # no. of pre-treatment periods
  brep      = 20,                  # no. of bootstraps (required for SEs)
  cluster   = 'country_id',        # variable to cluster SEs on
  parallel  = FALSE                 # run the bootstraps in parallel
)


tidy_m3 = tidy.did_multiplegt(m3)
tidy_m3


tidy_m3 |>
  within({
    term = gsub("^placebo_", "-", term)
    term = gsub("^effect", "0", term)
    term = gsub("^dynamic_", "", term)
    term = as.integer(term)
  }) |>
  ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange() +
  geom_hline(yintercept = 0, linetype = 'dashed') +
  labs(
    x = "Time to treatment (AF Growth Episode)", y = "Effect size", 
    title = "", 
    subtitle = ""
  )

ggsave("figures/Figure_3_AF_Growth_on_Democracy.pdf", height = 10, width = 20, units = "cm", dpi = 900)
ggsave("figures/Figure_3_AF_Growth_on_Democracy.png", height = 10, width = 20, units = "cm", dpi = 900)


############################################################################################
#### Figure F4: Analyzing Academic Freedom Growth's Effect on Democracy                ####
############################################################################################

m3b = did_multiplegt(
  vdem_subset_onset_sample, 
  Y = 'v2x_polyarchy', 'country_id', 'year', 'increase_episode', 
  controls = c("e_gdppc","e_pop", "v2x_jucon", "v2xlg_legcon"),
  dynamic   = 5,                  # no. of post-treatment periods
  placebo   = 10,                  # no. of pre-treatment periods
  brep      = 20,                  # no. of bootstraps (required for SEs)
  cluster   = 'country_id',        # variable to cluster SEs on
  parallel  = FALSE                 # run the bootstraps in parallel
)


tidy_m3b = tidy.did_multiplegt(m3b)
tidy_m3b


tidy_m3b |>
  within({
    term = gsub("^placebo_", "-", term)
    term = gsub("^effect", "0", term)
    term = gsub("^dynamic_", "", term)
    term = as.integer(term)
  }) |>
  ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange() +
  geom_hline(yintercept = 0, linetype = 'dashed') +
  labs(
    x = "Time to treatment (AF Growth Episode)", y = "Effect size", 
    title = "", 
    subtitle = ""
  )

ggsave("figures/Figure_3_AF_Growth_on_Democracy_with_covariates.pdf", height = 10, width = 20, units = "cm", dpi = 900)


############################################################################################
#### Figure 9: Analyzing Academic Freedom Decline Effect on Democracy                   ####
############################################################################################

m4 = did_multiplegt(
  vdem_subset_onset_sample, 
  Y = 'v2x_polyarchy', 'country_id', 'year', 'decline_episode', 
  dynamic   = 10,                  # no. of post-treatment periods
  placebo   = 10,                  # no. of pre-treatment periods
  brep      = 20,                  # no. of bootstraps (required for SEs)
  cluster   = 'country_id',        # variable to cluster SEs on
  parallel  = FALSE                 # run the bootstraps in parallel
)


tidy_m4 = tidy.did_multiplegt(m4)
tidy_m4


tidy_m4 |>
  within({
    term = gsub("^placebo_", "-", term)
    term = gsub("^effect", "0", term)
    term = gsub("^dynamic_", "", term)
    term = as.integer(term)
  }) |>
  ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange() +
  geom_hline(yintercept = 0, linetype = 'dashed') +
  labs(
    x = "Time to treatment (AF Decline Episode)", y = "Effect size", 
    title = "", 
    subtitle = ""
  )

ggsave("figures/Figure_4_AF_Decline_on_Demcoracy.pdf", height = 10, width = 20, units = "cm", dpi = 900)
ggsave("figures/Figure_4_AF_Decline_on_Demcoracy.png", height = 10, width = 20, units = "cm", dpi = 900)


############################################################################################
#### Figure F5: Analyzing Academic Freedom Decline Effect on Democracy                  ####
############################################################################################


m4b = did_multiplegt(
  vdem_subset_onset_sample, 
  Y = 'v2x_polyarchy', 'country_id', 'year', 'decline_episode', 
  controls = c("e_gdppc","e_pop", "v2x_jucon", "v2xlg_legcon"),
  dynamic   = 5,                  # no. of post-treatment periods
  placebo   = 10,                  # no. of pre-treatment periods
  brep      = 20,                  # no. of bootstraps (required for SEs)
  cluster   = 'country_id',        # variable to cluster SEs on
  parallel  = FALSE                 # run the bootstraps in parallel
)


tidy_m4b = tidy.did_multiplegt(m4b)
tidy_m4b


tidy_m4b |>
  within({
    term = gsub("^placebo_", "-", term)
    term = gsub("^effect", "0", term)
    term = gsub("^dynamic_", "", term)
    term = as.integer(term)
  }) |>
  ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange() +
  geom_hline(yintercept = 0, linetype = 'dashed') +
  labs(
    x = "Time to treatment (AF Decline Episode)", y = "Effect size", 
    title = "", 
    subtitle = ""
  )

ggsave("figures/Figure_4_AF_Decline_on_Demcoracy_with_covariates.pdf", height = 10, width = 20, units = "cm", dpi = 900)


############################################################################################
#### Supplementary Appendix: F6                                                         ####
#### Analyzing Autocratization Onset as a Function of Academic Freedom Decline Episodes ####
############################################################################################

vdem_subset_onset_sample <- vdem_subset_onset_sample %>%
  group_by(e_regionpol, year) %>%
  mutate(regionpol_EDI = mean(v2x_polyarchy, na.rm = TRUE)) %>%
  ungroup() %>%
  group_by(country_id) %>%
  mutate(e_gdppc_growth = ((e_gdppc - dplyr::lag(e_gdppc))/dplyr::lag(e_gdppc))*100) %>%
  mutate(e_gdppc_growth_roll5 = rollmean(e_gdppc_growth, k = 5, fill = NA)) %>%
  ungroup() 

## lags for independent variables ##

vdem_subset_onset_sample <- vdem_subset_onset_sample %>%
  mutate(e_gdppclog = log(e_gdppc+1), 
         e_poplog = log(e_pop +1), 
         year_1900 = year -1900) %>%
  group_by(country_id) %>%
  mutate(e_gdppclog_lag = lag(e_gdppclog), 
         e_gdppc_growth_roll5_lag = lag(e_gdppc_growth_roll5), 
         regionpol_EDI_lag = lag(regionpol_EDI), 
         e_poplog_lag = lag(e_poplog), 
         v2x_jucon_lag = lag(v2x_jucon), 
         v2xlg_legcon_lag = lag(v2xlg_legcon)) %>%
  mutate(region_factor = as.factor(e_regionpol_6C),
         region_factor = relevel(region_factor, ref=1), 
         year_sq = year_1900^2)


## drop NA ## 
vdem_subset_onset_sample <- vdem_subset_onset_sample %>%  
  group_by(country_id) %>%
  mutate(aut_ep_onset = ifelse(aut_ep==1 & lag(aut_ep)== 1, 0, aut_ep), 
         decline_episode_onset = ifelse(decline_episode==1 & lag(decline_episode)== 1, 0, decline_episode)) %>%
  mutate(decline_episode_onset_1 = lag(decline_episode_onset), 
         decline_episode_onset_2 = lag(decline_episode_onset, 2),
         decline_episode_onset_3 = lag(decline_episode_onset, 3),
         decline_episode_onset_4 = lag(decline_episode_onset, 4),
         decline_episode_onset_5 = lag(decline_episode_onset, 5),
         decline_episode_onset_1_3 = ifelse(decline_episode_onset_1 ==1 | decline_episode_onset_2 == 1 | 
                                              decline_episode_onset_3 == 1, 1, 0), 
         decline_episode_onset_1_5 = ifelse(decline_episode_onset_1 ==1 | decline_episode_onset_2 == 1 | 
                                              decline_episode_onset_3 == 1 | decline_episode_onset_4 == 1 |
                                              decline_episode_onset_5 == 1, 1, 0)) %>%
  dplyr::select(-c(decline_episode_onset_2, decline_episode_onset_3, decline_episode_onset_4, decline_episode_onset_5)) %>%
  ungroup() %>%
  filter(year >=1900) %>%
  drop_na(e_gdppclog, e_gdppc_growth_roll5_lag, regionpol_EDI_lag, 
          e_poplog,region_factor ,v2x_jucon_lag, v2xlg_legcon_lag, year_1900, country_id)  

#---------------------------------------------------#
#--------------------- Table F3    -----------------#
#---------------------------------------------------#

mod.1 <- brglm::brglm(as.factor(aut_ep_onset) ~  decline_episode_onset_1  + 
                        e_gdppclog_lag + e_gdppc_growth_roll5_lag +  e_poplog_lag +
                        regionpol_EDI_lag + v2x_jucon_lag + v2xlg_legcon_lag +
                        region_factor +
                        year_1900 + year_sq, 
                      data = vdem_subset_onset_sample, family = binomial(link = "probit"),
                      method = "brglm.fit", control.brglm = brglm.control(br.maxit = 1000, br.trace = T),
                      pl=T)

summary(mod.1)


mod.2 <- brglm::brglm(as.factor(aut_ep_onset) ~  decline_episode_onset_1_3  + 
                        e_gdppclog_lag + e_gdppc_growth_roll5_lag +  e_poplog_lag +
                        regionpol_EDI_lag + v2x_jucon_lag + v2xlg_legcon_lag +
                        region_factor +
                        year_1900 + year_sq, 
                      data = vdem_subset_onset_sample, family = binomial(link = "probit"),
                      method = "brglm.fit", control.brglm = brglm.control(br.maxit = 1000, br.trace = T),
                      pl=T)

summary(mod.2)



mod.3 <- brglm::brglm(as.factor(aut_ep_onset) ~  decline_episode_onset_1_5  + 
                        e_gdppclog_lag + e_gdppc_growth_roll5_lag +  e_poplog_lag +
                        regionpol_EDI_lag + v2x_jucon_lag + v2xlg_legcon_lag +
                        region_factor +
                        year_1900 + year_sq, 
                      data = vdem_subset_onset_sample, family = binomial(link = "probit"),
                      method = "brglm.fit", control.brglm = brglm.control(br.maxit = 1000, br.trace = T),
                      pl=T)

summary(mod.3)


modelsummary::modelsummary(list(mod.1, mod.2, mod.3),
                           output = "tables/Table1.tex", 
                           fmt = 3,
                           vcov = "stata", cluster = "country_id", 
                           estimate  = c("{estimate}{stars}"),
                           statistic = 'conf.int', 
                           conf_level = .95, 
                           coef_rename = c("(Intercept)" = "Intercept", 
                                           "decline_episode_onset_1" = "AF Decline Onset t-1", 
                                           "decline_episode_onset_1_3" = "AF Decline Onset (t-1 and t-3)", 
                                           "decline_episode_onset_1_5" = "AF Decline Onset (t-1 and t-5)", 
                                           "e_gdppclog_lag"  = "GDP pc log", 
                                           "e_gdppc_growth_roll5_lag" = "GDP growth", 
                                           "e_poplog_lag" = "Population log", 
                                           "regionpol_EDI_lag" = "Regional democracy levels", 
                                           "v2x_jucon_lag" =  "Judicial constraints on executive", 
                                           "v2xlg_legcon_lag" = "Legislative constraints on executive", 
                                           "region_factor5" = "Western Europe and North America", 
                                           "region_factor4" = "Subsaharan Africa", 
                                           "region_factor6" = "Asia and Pacific", 
                                           "region_factor1" = "Eastern Europe and Central Asia", 
                                           "region_factor3" = "MENA", 
                                           "year_1900" = "Year", 
                                           "year_sq" = "Year squared"))


modelsummary::modelsummary(list(mod.1, mod.2, mod.3),
                           fmt = 3,
                           vcov = "stata", cluster = "country_id", 
                           estimate  = c("{estimate}{stars}"),
                           statistic = 'conf.int', 
                           conf_level = .95, 
                           coef_rename = c("(Intercept)" = "Intercept", 
                                           "decline_episode_onset_1" = "AF Decline Onset t-1", 
                                           "decline_episode_onset_1_3" = "AF Decline Onset (t-1 and t-3)", 
                                           "decline_episode_onset_1_5" = "AF Decline Onset (t-1 and t-5)", 
                                           "e_gdppclog_lag"  = "GDP pc log", 
                                           "e_gdppc_growth_roll5_lag" = "GDP growth", 
                                           "e_poplog_lag" = "Population log", 
                                           "regionpol_EDI_lag" = "Regional democracy levels", 
                                           "v2x_jucon_lag" =  "Judicial constraints on executive", 
                                           "v2xlg_legcon_lag" = "Legislative constraints on executive", 
                                           "region_factor5" = "Western Europe and North America", 
                                           "region_factor4" = "Subsaharan Africa", 
                                           "region_factor6" = "Asia and Pacific", 
                                           "region_factor1" = "Eastern Europe and Central Asia", 
                                           "region_factor3" = "MENA", 
                                           "year_1900" = "Year", 
                                           "year_sq" = "Year squared"))


############################################################################################
#### Figure F6: Marginal Effects                                                        ####
############################################################################################

# Plot marginal effects 

library(ggeffects)
library(ggokabeito)
m1_pred <- ggpredict(mod.1, terms = c("decline_episode_onset_1[0,1]"),
                     type = "fe",
                     cov.fun = "vcovCL", 
                     vcov.type = "HC1",
                     vcov.args = list(cluster = vdem_subset_onset_sample$country_id)) %>%
  rename(af_decline = x)

m1_pred$af_decline <- as.factor(m1_pred$af_decline )

f1 <- ggplot(m1_pred, aes(x = af_decline, y = predicted, group = af_decline, color = af_decline, fill = af_decline)) +
  geom_point() +
  geom_linerange(aes(ymin = conf.low, ymax=  conf.high)) + 
  theme_bw() +
  ggtitle("") +
  labs(y= "P(Autocratization Onset)",
       x = "AF Decline Episode") + 
  theme(legend.position = "none") +
  scale_color_okabe_ito(order = c(2, 6)) +
  scale_fill_okabe_ito(order = c(2, 6)) +
  scale_y_continuous(limits = c(0, 0.25)) +
  scale_x_discrete(labels = c("0" = "FALSE", 
                              "1" = "TRUE"))

m2_pred <- ggpredict(mod.2, terms = c("decline_episode_onset_1_3[0,1]"),
                     type = "fe",
                     cov.fun = "vcovCL", 
                     vcov.type = "HC1",
                     vcov.args = list(cluster = vdem_subset_onset_sample$country_id)) %>%
  rename(af_decline = x)

m2_pred$af_decline <- as.factor(m2_pred$af_decline )

f2 <- ggplot(m2_pred, aes(x = af_decline, y = predicted, group = af_decline, color = af_decline, fill = af_decline)) +
  geom_point() +
  geom_linerange(aes(ymin = conf.low, ymax=  conf.high)) + 
  theme_bw() +
  ggtitle("") +
  labs(y= "P(Autocratization Onset)",
       x = "AF Decline Episode") + 
  theme(legend.position = "none") +
  scale_color_okabe_ito(order = c(2, 6)) +
  scale_fill_okabe_ito(order = c(2, 6)) +
  scale_y_continuous(limits = c(0, 0.25)) +
  scale_x_discrete(labels = c("0" = "FALSE", 
                              "1" = "TRUE"))


m3_pred <- ggpredict(mod.3, terms = c("decline_episode_onset_1_5[0,1]"),
                     type = "fe",
                     cov.fun = "vcovCL", 
                     vcov.type = "HC1",
                     vcov.args = list(cluster = vdem_subset_onset_sample$country_id)) %>%
  rename(af_decline = x)

m3_pred$af_decline <- as.factor(m3_pred$af_decline )

f3 <- ggplot(m3_pred, aes(x = af_decline, y = predicted, group = af_decline, color = af_decline, fill = af_decline)) +
  geom_point() +
  geom_linerange(aes(ymin = conf.low, ymax=  conf.high)) + 
  theme_bw() +
  ggtitle("") +
  labs(y= "P(Autocratization Onset)",
       x = "AF Decline Episode") + 
  theme(legend.position = "none") +
  scale_color_okabe_ito(order = c(2, 6)) +
  scale_fill_okabe_ito(order = c(2, 6)) +
  scale_y_continuous(limits = c(0, 0.25)) +
  scale_x_discrete(labels = c("0" = "FALSE", 
                              "1" = "TRUE"))


combined_plot <- ggarrange(f1, f2, f3, ncol = 3, labels= c("1 year lag", "1-3 year lag", 
                                                           " 1-5 year lag"))

ggsave("figures/Autocratization_Onset_Marginal.pdf", combined_plot, width = 18, height = 10, units = "cm", dpi = 900, device = "pdf")

